% Model function file for social planner's optimal harp seal harvest problem
% s(:,1) -> state variable (harp seal stock, S)
% s(:,2) -> state variable (mean sea ice, Z)
% x -> control variable (harvest, q)
% Vs(:,1) -> partial derivative of value function wrt resource stock S
% Vs(:,2) -> partial derivative of value function wrt sea ice Z
function [out1,out2]=mff2varharp_final(flag,s,x,Vs,Vss,p0,p1,c0,c1,u0,delta,lambda,zeta,K1,b,m0,m1,gamma,sig0,sig1,scale,scale2,theta,D,Zmax)
switch flag
case 'x' % optimal value of the control variable given S and Vs
  Cost=c0*scale-c1*scale*s(:,2); % marginal cost function (%/seal harvested)
  out1=max((p0*scale-(Cost+scale2*Vs(:,1)-D*0.5*scale2*Vss(:,1)/scale))./(p1*s(:,1)*scale^2),0);  
case 'f' % objective function given S and x
    Cost=c0*scale-c1*scale*s(:,2); % marginal cost function
    out1=(x.*s(:,1)*scale*p0-0.5*p1*(x.*s(:,1)*scale).^2-Cost.*x.*s(:,1)+u0*(log(scale)+log(s(:,1))))/scale2; % social surplus function
case 'g' % case 'g' desingates the drift rate of the stochastic process
    Growth=s(:,1).*b.*(1-m0+m1*s(:,2)).*(1-lambda*scale*s(:,1))-zeta*s(:,1);  
    out1=[Growth-x.*s(:,1), -gamma*s(:,2)];
case 'sigma'
  out1(:,1,1)=(s(:,1).^2.*b^2.*m1^2.*(sig0-(theta*sig1*s(:,2)+(1-theta)*sig1*Zmax)).*(1-lambda*scale*s(:,1)).^2+D*(s(:,1).*b/scale.*(1-m0+m1*s(:,2)).*(1-lambda*scale*s(:,1))+zeta/scale*s(:,1)+x./scale.*s(:,1))).^0.5;
  %out1(:,1,1)=0.1*s(:,1);
  out1(:,2,2)=0*s(:,2);
  out1(:,2,1)=0;
case 'rho'
  out1=delta+zeros(size(s(:,1),1),1);  
end